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1.  SUMMARY 


The  purpose  of  this  project  was  to  develop  both  linear  and  nonlinear  numerical  modeling 
codes  for  infrasound  propagation  and  to  use  them  to  compute  propagation  through 
realistic  atmospheric  models.  The  objective  was  to  determine  whether  nonlinear  effects 
have  a  bearing  on  estimates  of  source  magnitudes  derived  using  either  infrasound 
amplitudes  or  dominant  frequencies.  Nonlinearity  is  expected  to  primarily  affect  the  near¬ 
source  regions  for  large  explosions  or  naturally  occurring  events,  and  inffasound 
propagation  in  the  thermosphere,  where  ambient  pressures  are  extremely  low.  Other  areas 
in  which  nonlinearity  may  play  a  role  is  in  the  infilling  of  acoustic  shadow  zones,  and 
propagation  near  caustics. 

The  Navier-Stokes  equations  are  solved  using  a  finite-difference,  time-domain  (FDTD) 
approach  for  axi-symmetric  environmental  models,  allowing  three-dimensional  acoustic 
propagation  to  be  simulated  using  a  two  dimensional  Cylindrical  coordinate  system.  A 
method  to  stabilize  the  FDTD  algorithm  in  a  viscous  medium  at  atmospheric  densities 
characteristic  of  the  lower  thermosphere  is  described.  The  stabilization  scheme  slightly 
alters  the  governing  equations,  but  results  in  quantifiable  dispersion  characteristics.  It  is 
shown  that  this  method  leaves  sound  speeds  and  attenuation  unchanged  at  frequencies 
that  are  well  resolved  by  the  temporal  sampling  rate,  but  strongly  attenuates  higher 
frequencies.  Numerical  experiments  are  performed  to  assess  the  effect  of  source  strength 
on  the  amplitudes  and  spectral  content  of  signals  recorded  at  ground  level  at  a  range  of 
distances  from  the  source.  It  is  shown  that  the  source  amplitudes  have  a  stronger  effect  on 
a  signal’s  dominant  frequency  than  on  its  amplitude.  Applying  the  stabilized  code  to 
infrasound  propagation  through  realistic  atmospheric  profiles  shows  that  nonlinear 
propagation  alters  the  spectral  content  of  low  amplitude  thermospheric  signals, 
demonstrating  that  nonlinear  effects  are  significant  for  all  detectable  thermospheric 
returns. 

2.  INTRODUCTION 

One  of  the  primary  drivers  of  infrasound  research  in  the  past  two  decades  has  been  its 
application  to  nuclear  monitoring  (Brachet  et.al.,  2009).  A  key  goal  of  this  research  has 
been  to  find  a  reliable  method  to  estimate  the  source  yield  of  large  explosions  based  on 
infrasound  recordings.  Empirical  source  yield  relations  have  been  derived  for  the 
dominant  period  of  an  infrasound  signal  (Revelle,  1997;  Ens,  2012)  and  for  its  amplitude 
(Blanc  et.al,  1997;  ANSI,  1983)  based  on  observations  of  a  large  number  of  nuclear 
explosions  or  bolides.  Large  variations  in  the  yield  estimates  remain  because  these 
relations  do  not  take  into  account  sound  and  wind  speed  variability,  which  strongly 
affects  infrasound  propagation  (e.g.  de  Groot-Hedlin  et.al.,  2010).  Le  Pichon  et.al.  (2012) 
used  linear  parabolic  equation  (PE)  methods  to  provide  more  realistic  expressions  for 
inffasound  transmission  loss  (TL)  for  a  wide  range  of  realistic  atmospheric  sound  speed 
profiles.  These  expressions  provide  an  improved  description  of  the  physics  of  long-range 
propagation  and  are  valid  for  both  stratospheric  and  thermospheric  ducting. 
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A  better  description  of  infrasound  propagation,  especially  for  large  explosive  sources, 
requires  incorporating  nonlinear  effects  into  numerical  modeling  algorithms.  Nonlinear 
propagation  implies  that  the  amplitude  of  the  pressure  and  density  perturbations  due  to 
the  passage  of  a  sound  wave  is  a  non-negligible  fraction  of  the  ambient  fields. 
Nonlinearity  affects  both  the  amplitude  and  frequency  of  infrasound  signals  and  therefore 
impacts  source  yield  estimates.  Since  nonlinear  effects  are  dependent  on  the  absolute 
amplitude  of  acoustic  arrivals,  and  because  acoustic  amplitudes  decrease  much  more 
slowly  for  a  line  source  than  for  a  spherical  source,  numerical  algorithms  that  compute 
propagation  from  a  spherical  source  must  be  used  to  compare  linear  vs.  nonlinear 
infrasound  propagation.  The  finite  difference  time  domain  (FDTD)  method  developed 
here  is  similar  to  that  of  de  Groot-Hedlin  (2012)  in  that  the  sources  and  models  are 
required  to  be  axisymmetric.  This  requirement  allows  for  the  efficient  solution  of  the 
governing  equations  in  a  2D  cylindrical  coordinate  system,  with  quantities  that  may  vary 
in  both  the  radial  and  vertical  directions.  This  leads  to  a  FDTD  code  that  can  be  used  to 
simulate  3D  infrasound  propagation  from  a  spherical  source. 

In  this  paper,  both  linear  and  nonlinear  FDTD  propagation  modeling  codes  are  applied  to 
several  environmental  models  to  find  where  nonlinearity  affects  the  synthesized 
wavefields.  Nonlinearity  affects  infrasound  not  only  near  the  shock  front,  which  is  not 
examined  here,  but  also  affects  infrasound  propagation  within  the  thermosphere,  where 
ambient  pressures  are  extremely  low.  Other  areas  in  which  nonlinearity  may  play  a  role  is 
in  the  infilling  of  acoustic  shadow  zones,  and  propagation  near  caustics.  Ultimately,  the 
goal  is  to  determine  whether  nonlinear  effects  have  a  bearing  on  estimates  of  source 
magnitudes  derived  using  either  infrasound  amplitudes  or  dominant  frequencies. 
However,  a  straight-forward  application  of  the  FDTD  governing  equations  can  be 
unstable  in  the  presence  of  shock  fronts  or  in  viscous  media  at  densities  characteristic  of 
the  thermosphere. 

Several  approaches  have  been  developed  to  stabilize  FDTD  algorithms  in  the  presence  of 
discontinuities.  Walkington  and  Eversman  (1986)  introduced  an  artificial  viscosity  term 
to  dissipate  oscillations  that  are  poorly  resolved  by  the  spatial  grid,  resulting  in  a 
dissipation  term  that  varies  with  the  fourth  power  of  the  frequency,  in  addition  to  the  real 
viscous  attenuation  term  that  varies  with  the  square  of  the  frequency.  Sparrow  and  Raspet 
(1991)  use  this  technique  to  stabilize  simulations  of  finite  amplitude  acoustic 
propagation.  Marsden  et.al.  (2014)  employ  a  shock-capturing  filter  technique  that 
suppresses  Gibbs’  oscillations  that  arise  near  shock  waves.  It  operates  by  applying  a 
nonlinear  low  pass  filter  at  each  time  step  (Bogey  et.al.,  2009).  Shepherd  et.al.  (2009)  use 
the  weighted  essentially  non-oscillatory  (WENO)  method  to  stably  propagate  acoustic 
shocks.  The  WENO  scheme  is  another  shock-capturing  method  that  uses  adaptive  stencils 
to  avoid  interpolation  across  discontinuities  (Shu,  1998)  to  avoid  introducing  oscillations 
into  the  solution.  However,  it  can  have  the  unintended  effect  of  excessively  damping 
relevant  small-scale  structure  (Bogey  et.al,  2009).  Moreover,  the  WENO  method  is  very 
computationally  intensive  as  several  stencils  are  computed  at  each  iteration.  Costa  and 
Don  (2007)  employ  a  hybrid  approach  that  uses  the  WENO  scheme  near  discontinuities, 
and  a  central  finite  difference  method,  coupled  with  high  order  filtering,  that  strongly 
dissipates  high  frequency  oscillations  over  the  remainder  of  the  grid. 
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Most  FDTD  stabilization  method  work  by  adding  steps  to  the  basic  differencing  schemes 
and  in  the  process  they  alter  the  governing  equations  being  modeled.  This  has  the  desired 
effect  of  attenuating  high  frequency  fluctuations  that  are  poorly  resolved  by  the  spatial 
grid  but  also  alters  the  propagation  velocity,  although  this  side  effect  is  not  generally 
recognized.  These  changes  result  from  the  fact  that  modifying  the  governing  equations 
alters  the  wavenumber,  which  is  related  to  both  the  attenuation  and  sound  speed.  In  many 
stabilization  schemes,  the  changes  to  the  wavenumbers  are  not  easily  quantifiable.  In  this 
report,  a  new  method  to  improve  the  stability  of  the  FDTD  algorithm  is  presented.  Aside 
from  being  significantly  less  computationally  intensive  than  other  stabilization  methods, 
an  advantage  is  that  the  frequency  dependent  effect  on  the  sound  speed  profile  and 
attenuation  can  be  quantified. 

This  report  is  organized  as  follows.  The  sets  of  equations  governing  both  linear  and 
nonlinear  acoustic  wave  propagation  are  derived  in  the  Sec  3.1.  In  section  3.2,  a  robust 
FDTD  approach  to  solving  these  equations  is  presented,  along  with  a  dispersion  analysis 
of  the  effect  of  the  stabilization  terms  on  the  governing  equations.  Several  numerical 
examples  are  presented  in  Sec.  4  to  show  how  the  spectral  content  and  amplitudes  of 
infrasound  signals  are  altered  as  they  propagate  away  from  the  source,  and  the  effect  of 
increasing  source  amplitude  on  these  characteristics.  The  report  concludes  with  remarks 
on  the  implications  of  the  results  for  nuclear  monitoring  and  to  furthering  basic 
understanding  of  atmospheric  properties  within  the  thermosphere. 


3.  TECHNICAL  APPROACH 


3.1  Governing  Equations  for  Infrasound  Propagation 


The  fluid  equations  governing  infrasound  propagation  are  derived  from  the  conservation 
of  mass, 


Dp 

Dt 


-pV»v, 


(1) 


the  conservation  of  momentum  which  is  described  by  the  Navier  Stokes  equation  for  a 
viscous  fluid, 


p—  =  -Vp+pV2v +{pB  +  p/  3)V(V  •  v), 


and  the  conservation  of  energy 


,V-rtV.v)+ 


dv  dv  dv  dv 

T  -  +  T  -  +  T  ~+T  - 

"dr  "  dz  zr  dr  22  dz 


(2) 


(3) 


where  the  internal  energy  is  given  by 
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(4) 


e  = 


This  is  an  extended  set  of  Navier-Stokes  equations  that  relates  the  time  rate  of  change  of 
the  density  p  and  particle  velocity  vand  the  internal  energy  e,  to  the  total  pressure,  p,  and 

the  shear  and  bulk  viscosities  fj,  and  fiR.  The  relevant  components  of  the  stress  tensor  are 
given  by 


T  =fl 


dv 

> _ r 

dr 


T  =T  =U 

1-7.  7.r  * 

dv 


t  =  a 

77  4 


dz 


“I™. 

(5a) 

dv  dv 

_ L-| _ £. 

dz  dr 

(5b) 

-fcv*p) 

(5c) 

7 

Dt’  n 

In  Eqn.  (2),  only  the  irrotational  part  of  the  particle  velocity  is  relevant  because  the 
density  is  related  to  the  divergence  of  the  particle  velocity,  so  it  may  be  re-written  as 


Dv  ~ 

/?— =  -V/;+//,V  v, 

Dt 


(6) 


where  //  is  the  total  viscosity  fig  +  //  /  3 .  In  this  study,  the  value  for  /i  in  air  is  set  to  1.8 

x  10'5  kg  m'V1  (Pierce,  1981)  and  the  bulk  viscosity  is  set  to  zero.  More  realistically,  the 
dynamic  viscosity  is  a  function  of  temperature  (Sutherland  and  Bass,  2004);  it  varies  by 
approximately  a  factor  of  two  for  realistic  atmospheric  temperatures  up  to  120  km. 
However,  given  uncertainties  in  attenuation  mechanisms  at  high  altitude  (Akintunde  and 
Petculescu,  2014),  the  viscosity  is  represented  by  a  constant  in  this  study. 

Equations  (1),  (3),  (4)  and  (6)  form  a  complete  set  of  governing  equations  for  acoustic 
propagation  in  a  viscous  medium.  Thermal  conductivity  and  molecular  relaxation  are 
neglected;  these  processes  have  little  effect  at  infrasound  frequencies  (Sutherland  and 
Bass,  2004).  The  equations  do  not  include  gravity  terms,  so  they  apply  to  frequencies 
above  the  natural  buoyancy  frequency  of  the  atmosphere  (the  Brunt- Vaisala  frequency), 
which  ranges  from  about  about  1  mHz  in  the  troposphere  to  3  mHz  in  the  stratosphere 
and  thermosphere  (e.g.  Lingevitch  et  al  1999)..  This  set  of  equations  is  similar  to  that 
developed  in  Wochner  et.  al.  (2005),  Shepherd  et.  al.  (2009),  and  de  Groot-Hedlin  (2012) 
except  that  the  energy  equation  is  used  instead  of  the  entropy  balance  equation. 

Equations  (1),  (3)  and  (4)  may  be  combined  to  show  that  the  pressure  and  density  are 
related  as 
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(7) 


—  =  c 


Dp  _  2Dp 
Dt~C  Dt' 


Dt 


for  linear  propagation  in  an  ideal  gas,  where  c  is  the  local  sound  speed.  Combining  Eqs. 
(1)  and  (7)  yields 


(8) 


Equations  (6)  and  (8)  form  a  complete  set  of  governing  equations  for  linear  acoustic 
propagation. 

In  the  absence  of  shear  viscosity  and  other  absorption  effect,  i.e.  in  an  isentropic  process 
in  an  ideal  gas,  eqns.  (1),  (3),  and  (4)  may  be  combined  to  show  that  the  pressure  is 
related  to  the  density  as 


(9) 


where  the  o  subscripts  indicate  ambient  values,  and  jds  the  specific  heat  ratio,  which  is 
1.4  in  air.  Equations  (1),  (6)  and  (9)  form  a  complete  set  of  governing  equations  for 
acoustic  propagation  in  an  isentropic  medium. 

In  what  follows,  the  second  set  of  equations  (6  and  8)  is  used  to  derive  the  coupled 
differential  equations  for  linear  acoustic  propagation,  and  the  third  set  (Eqs.  (1),  (6)  and 
(9))  is  used  to  derived  the  equivalent  expressions  for  nonlinear  propagation  for  a  medium 
without  viscosity.  The  advantage  of  using  the  latter  set  of  equations  for  the  nonlinear 
problem  is  that  we  avoid  direct  computation  of  the  perturbed  sound  speed.  For  the  linear 
problem,  perturbations  in  the  sound  speeds  due  to  the  passage  of  an  acoustic  pulse  are 
included  in  second  and  higher  order  terms  and  are  neglected.  Finally,  the  first  set  of 
equations  (1,  3,  4,  and  6)  is  used  to  compute  nonlinear  propagation  through  a  viscous 
atmosphere. 

The  following  derivations  for  the  equations  governing  both  linear  and  nonlinear 
propagation  make  the  assumption  that  the  ambient  temperature,  pressure  and  density 
fields  vary  only  with  altitude.  This  is  a  realistic  representation  of  the  atmosphere  to  first 
order,  as  ambient  density  and  pressure  variations  are  due  to  the  effects  of  gravity.  Since 
the  ambient  densities  and  pressures  are  stratified  by  gravity  even  though  gravity  terms  are 
not  explicitly  included  in  this  formulation,  terms  involving  vertical  gradients  of  the 
ambient  pressure  and  density  fields  are  omitted.  Wind  is  also  neglected  in  this  study  in 
order  to  more  accurately  isolate  the  effects  of  differences  in  linear  vs.  nonlinear 
propagation.  An  axisymmetric  cylindrical  coordinate  system  is  used  for  each 
formulation,  leading  to  a  code  that  may  be  used  to  compute  inffasound  propagation  from 
a  spherical  source  embedded  in  a  realistic  atmosphere. 
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A.  Linear  propagation: 

If  the  ambient  pressure  p0,  density,  p0,  and  temperature  vary  only  with  altitude  and  are  in 
hydrostatic  equilibrium,  then  Eqs.  (8)  and  (6)  may  be  linearized  as 


dp  J  dv  dv  v 

-L±  =  -p(?\  — r-+— z-+-^\, 

00 1,  dr  dz  rj 


and 


dt 


dv 

_ r_ 

dt 

dv 

_ z_ 

dt 


1  fy  +  n 


p  dr  p 

*  n  *  t 


o\ 


d2v  1  dv  d2v  v  ^ 

- T  + - -  + - 
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1  SP,  M, 


d2v  1  dv  d2v 
"+ - -  +  - 


pn  dz  po  dr2  r  dr  dz2  j 


(10) 


(11a) 


(lib) 


where  p  denotes  the  perturbed  pressure  and  the  perturbed  particle  velocity  has  been 
written  asv  =  [v„0,vz].  This  is  a  valid  approximation  as  long  as  the  perturbed  quantities 
are  a  negligible  fraction  of  the  ambient  quantities. 


B.  Nonlinear  propagation  in  an  inviscid  medium: 

In  the  nonlinear,  inviscid  case,  perturbed  quantities  may  be  large  so  higher  order  terms 
are  retained.  In  this  case,  Eq.  (1)  becomes 


dJ\ 
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=  ~{P+P, 


f  dv  v  dv  ^  ^ 
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and  Eq.  (6)  becomes 


The  equation  for  pressure  perturbations  is  derived  from  Eq.  (9)  as 

V  pY 

P'=P"  l1+*J  -1 


(12) 
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—  v  — -+v  — -  , 
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(13b) 

(14) 


Approved  for  public  release;  distribution  is  unlimited. 
6 


C.  Nonlinear  propagation  in  a  viscous  medium: 

For  the  more  general  case  of  nonlinear  propagation  in  a  viscous  medium,  the  momentum 
equations  are 


dvr  1 

dPs 

f  r\2 

a  v 

r 

dt  P  +P 

ar  ^ 

Kdr2 

dv  1 

z  _ 

1 

1 

fcfv 

z 

dt  p  +p 

1  o  1  s 

dz  c 

l dr 2 

vY 


1  dv  d2v 

- -  + — f- — I  -I  v 

r  dr  dz2  r  )  {  ' 


( 


dv 

— — +v 
r  dr  2  dz 


dv^ 1 

f  J'(15a) 


1  dv  dV  ^ 

■  + - ~+ — f 

r  dr  dz  j 


f  dv  dv  1 


lv' 


r  dr 


- + v 


dz 


)■ 
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The  equation  for  internal  energy,  Eqn.  3  is  given  by 

de  p  +p  ( dv  v  dv  \  (  de  de  ^ 

-  =  — 5  -  -+—+ — -  -  v  -  +  v  -\+A, 
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(16) 


where  the  attenuation  term,  A,  is  given  by 
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Finally,  the  pressure  perturbation  is  derived  from  the  Eqn.  4  as 


p  ={e  +e  )(y-l)(/>  +p  )-p  . 

5  05  0.0  (18) 

The  conservation  of  mass  is  as  given  by  Eq.  (12)  .This  formulation  could  also  be  used 
for  nonlinear  acoustic  propagation  for  an  inviscid  medium,  however  the  equations  in 
section  3.1  are  more  computationally  efficient.  Similarly,  the  equations  in  this  section 
could  be  linearized  for  use  in  linear  propagation  modeling  but  the  formulation  in  3. 1  is 
more  efficient. 


3.2  The  FDTD  Numerical  Algorithm 

A.  The  FDTD  Grids 

Each  set  of  equations  in  the  previous  section  is  a  complete  set,  in  a  form  amenable  to 
solution  by  the  FDTD  methods.  Finite  difference  methods  yield  solutions  to  differential 
equations  by  replacing  derivatives  of  continuous  functions  by  their  finite  difference 
approximations  formed  over  sets  of  discrete  nodes. 

A  spatially  staggered  grid,  as  shown  in  Figure  1  is  used  to  implement  the  FDTD 
algorithms.  The  solution  domain  is  divided  into  a  grid  of  square  cells  of  dimension  A, 
each  with  a  uniform  density  p° ,  and  static  sound  speed  c,  with  the  axisymmetric  source 
placed  at  r=0  along  the  left  boundary.  Pressure  variables  are  sampled  at  each  cell  center 
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and  particle  velocities  are  sampled  along  the  cell  boundaries.  The  primary  advantage  of 
using  the  spatially  staggered  grid  is  that  it  provides  second  order  accurate  central 
differences  over  a  single  node  for  spatial  derivatives  in  the  FDTD  implementation. 
Higher  order  spatial  differences  could  lead  to  equivalent  accuracy  with  a  coarser  grid 
sampling  in  the  smooth  regions  of  the  solution.  However,  this  leads  to  a  loss  in  accuracy 
where  the  solution  changes  rapidly,  which  can  result  in  instability.  The  problem  of 
numerical  stability  is  discussed  further  in  section  3.2B. 
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Figure  1.  The  FDTD  model  grid  setup.  The  spatial  model  is  discretized  as  a  set  of  cells, 
indicated  by  the  solid  lines,  each  with  uniform  ambient  density p ,  and  static  sound  speed 
c.  Pressure  and  density  perturbations  are  defined  at  nodes  at  each  cell  center,  shown  by 
circles.  The  horizontal  particle  velocity  vr  (triangles  pointing  right)  and  vertical  particle 
velocity  variables  vz  (triangles  pointing  up)  nodes  straddle  the  cell  boundaries  of  the 
spatially  staggered  grid.  The  solution  domain,  marked  by  the  thick  line,  is  bounded  by 
a  rigid  surface  at  z=0,  and  by  an  axis  of  symmetry  at  r=0.  PML  boundaries  are  used 
to  simulate  open  boundaries  at  the  top  and  right  edges  of  the  model. 


The  spatial  solution  domain  is  bounded  at  the  bottom  by  a  rigid  surface  at  z=0,  and  at  left 
by  an  axis  of  symmetry  at  r=0.  The  rigid  surface  is  specified  by  setting  vz  (r,z=0)  =  0. 

The  vertical  component  of  particle  velocity  vz  is  anti-symmetric  and  other  terms  are 
symmetric  about  z=0.  Since  the  model  is  axially  symmetric,  the  radial  component  of 
particle  velocity  is  anti-symmetric  about  r=0  and  all  other  terms  are  symmetric.  This 
means  that  vr  (0,z)  =  0  at  the  left  edge  of  the  model,  a  result  that  can  be  confirmed  by 
applying  l’Hopital’s  rule  to  Eqs.  (11a) .  Several  terms  in  the  FDTD  approximation  to  the 
governing  equations  require  field  values  at  r<0  or  z<0,  beyond  the  model  boundaries. 
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This  is  resolved  by  making  use  of  field  symmetries.  For  instance,  the  first  order  spatial 
derivatives  that  appear  in  the  attenuation  terms  are  formed  over  adjacent  particle  velocity 

nodes.  Thus,  the  dv?  /drterm  at  r=A/2  requires  values  at  i/(r  =  -A/2,z)  and 

vz(r  =  3A/2,z)  to  form  a  central  difference  derivative.  By  symmetry,  vz  (-r,z)  =  vz  (r,z) 

so  that  the  dvzj dr  term  at  the  leftmost  nodes  is  computed  as 

dvz(A/2,Zj)  v.(3A/2,z,)-v?(-A  12, Zj)  v.pAI2,Zj)-vz(M2,Zj) 

- = - = - .  (19) 

dr  2  A  2A 

Where  particle  velocities  are  required  at  a  spatially  staggered  node,  the  particle  velocities 
are  averaged.  That  is,  where  vzis  required  at  a  vr  node  as  in  Eq.  (13a),  or  vr  at  a  vz  node 
as  in  Eq.  (13b),  the  four  nearest  particle  velocity  values  are  averaged. 

The  inclusion  of  attenuation  terms  in  the  linear  equations,  and  self-advection  term  in  the 
nonlinear  equations  complicate  the  standard  method  of  staggering  the  velocities  and 
pressures  in  both  space  and  time.  These  terms  require  that  spatial  and  time  derivatives  are 
computed  simultaneously,  rather  than  serially  as  in  the  staggered  leap-frog  approach 
(Yee,  1966).  Instead  a  “non-staggered  leap-frog”  scheme,  in  which  centered  temporal 
differences  are  formed  over  two  time  steps,  is  applied.  This  provides  second  order 
accuracy  in  the  time  derivatives  and  avoids  the  use  of  one-sided  differences  which  can 
introduce  attenuation  or  growth  factors  into  the  solution  (de  Groot-Hedlin,  2008).  This 
approach  was  previously  applied  by  Ostashev  et  al  (2005)  for  computing  acoustic  waves 
advected  through  a  windy  atmosphere,  by  de  Groot-Hedlin  (2008)  for  computing  acoustic 
propagation  through  a  viscous  atmosphere,  and  de  Groot-Hedlin  et.al.  (2011)  for 
computing  infrasound  propagation  through  a  windy,  viscous  atmosphere. 

The  grid  and  time  sampling  rates  in  an  FDTD  simulation  depend  on  the  highest 
frequencies  used  for  propagation  modeling.  Following  Taflove  and  Hagness  (2000),  the 
spatial  discretization  used  is  10  nodes  per  wavelength  at  the  lowest  sound  speed  and 
highest  frequency  of  interest  in  the  computations.  For  three  dimensional  acoustic 
propagation  in  an  inviscid  medium,  the  Courant  stability  limit  for  a  Yee-type  staggered 
leap-frog  approach  is  given  by  A t  <  A/(cm.lx  V3).  (Taflove  and  Hagness,  2000).  In  the 
algorithms  developed  here,  the  time  derivatives  are  formed  over  two  time  steps  so  each 
time  step  is  half  the  Courant  limit.  Because  energy  shifts  from  lower  to  higher 
frequencies  for  nonlinear  computations,  the  requirement  of  10  nodes  per  wavelength 
implies  finer  spatial  and  temporal  discretizations  are  needed  as  the  source  amplitude 
increases. 

A  perfectly  matched  layer  (PML)  absorbing  boundary  condition  (Berenger,  1994)  that 
strongly  absorbs  outgoing  waves  from  the  computational  region  is  imposed  on  the  top 
and  right  side  of  the  model.  Berenger’ s  original  PML  scheme  applies  to  a  Cartesian 
coordinate  system,  but  was  modified  for  use  in  cylindrical  coordinates  by  retaining  the 
split-field  formulation  (the  pressure  in  the  boundary  layer  is  split  into  vertically  and 
radially  propagating  components)  and  adding  cylindrical  terms  to  the  radial  part  of  the 
pressure  field  in  the  boundary  region.  The  PML  method  has  not  been  modified  to  account 
for  attenuation  or  nonlinearity. 
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B.  Numerical  stability  and  dispersion  relations 

Any  solution  derived  using  an  FDTD  method  is  inexact,  in  part  because  the  continuous 
differential  equations  are  replaced  by  discretized  approximations.  These  errors  can  be 
quantified  using  dispersion  relations  that  are  derived  by  considering  harmonic  travelling 
wave  solutions,  as  in  Taflove  and  Hagness  (2000)  for  electromagnetic  problems,  or  de 
Groot-Hedlin  et.al  (2011)  for  infrasound  propagation.  A  second  reason  is  that  the  finite 
precision  of  any  computer  allows  very  small  truncation  errors  to  creep  into  a  solution; 
over  a  very  large  number  of  time  steps  these  errors  can  accumulate,  causing  the  solution 
to  become  unstable.  A  third  source  of  error  is  the  presence  of  physical  discontinuities  in 
the  model  or  shock  waves  in  the  solution.  The  FDTD  method  works  by  discretizing  field 
variables  over  discrete  stencils,  or  sets  of  nodes.  An  increase  in  the  number  of  nodes  used 
to  approximate  the  derivative  leads  to  higher  orders  of  accuracy  and  allows  for  coarser 
grid  discretization,  so  long  as  the  function  is  smooth  within  the  width  of  the  stencil. 
However,  the  use  of  a  wide  stencil  in  regions  where  the  solution  undergoes  rapid  change  - 
such  as  near  shock  waves  or  physical  discontinuities  in  the  model  -  seriously  degrades  the 
accuracy,  and  produces  large  oscillations  in  the  solution.  These  oscillations  do  not 
decrease  with  increasing  grid  refinement  and  create  instabilities  in  the  FDTD  method.  In 
this  section,  a  method  is  presented  that  suppresses  instabilities  that  arise  in  regions  where 
the  kinematic  viscosity  is  high. 

A  straightforward,  unstabilized  implementation  of  the  FDTD  method  to  the  linear  set  of 
equations  involves  solving  for  the  acoustic  pressure  perturbation  as, 

pN+1  =pN  1 +  FpA(vN)2At,  (20) 

where  the  superscript  represents  the  number  of  time  steps,  and  the  notation  FPA  denotes 
the  finite  difference  operator  that  replaces  the  acoustic  term  in  Eq.  (10).  Equation  (1 1)  for 
the  particle  velocity  can  be  represented  as 

vN+1  =  1  +(F>")+  FJvN)>  2M  (21) 

where  FVA  and  FVfi  denote  the  finite  difference  operators  for  the  acoustic  term  and  the 
viscosity  terms,  respectively.  Expressions  for  the  nonlinear  density  and  particle 
velocities  can  be  represented  similarly.  These  equations  are  stable  for  inviscid  models, 
provided  that  the  spatial  and  temporal  discretizations  are  sufficiently  fine  to  accurately 
sample  the  acoustic  perturbations.  However,  the  addition  of  the  viscosity  term  renders 
these  equations  unstable  at  densities  consistent  with  those  in  the  thermosphere,  even  for 
linear  propagation  in  a  medium  with  realistic  viscosities.  This  instability  is  likely  due  to  a 
combination  of  computer  truncation  errors  and  the  fact  that  the  introduction  of  high 
kinematic  viscosities  can  decrease  the  time  step  required  to  maintain  stability.  A  full 
Fourier  analysis  that  takes  into  account  both  viscosity  and  acoustic  propagation  may  be 
required  to  find  criterion  for  a  stable  time  step.  Such  an  analysis  is  beyond  the  scope  of 
this  paper. 

Solving  for  the  particle  velocity  in  a  two-step  procedure  stabilizes  the  equations.  First, 
the  acoustic  velocity  is  added  to  provide  an  interim  field, 
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(22) 


v  =  vN-'  +  FVA(pN)2At 

then  the  particle  velocity  is  computed  by  applying  the  viscosity  operator  to  this  interim 
field: 


— * L  N  + 1  *  77 

y  =  v  +  Fx 


Vu 


(23) 


Combining  Eqs.  (22)  and  (23)  yields  a  stable  form  for  computing  particle  velocity: 


vN"  =  vN  1  +  FVA[pN)2At+  FVfi(vNl  +  FVA{pN)2At)2At .  (24) 

This  may  be  compared  to  the  unstabilized  form,  Eq.  (21).  The  algorithm  to  stabilize  the 
nonlinear  equations  is  similar;  the  acoustic  and  self-advection  terms  are  first  computed  to 
find  the  particle  velocity  for  the  inviscid  case,  then  the  attenuation  term  operates  on  the 
resulting  field.  This  method  allows  for  the  direct  computation  of  the  dispersion  relation 
connected  to  the  stabilized  form  of  the  FDTD  equations,  which  leads  to  the  effective 
sound  speed  and  attenuation  as  a  function  of  altitude.  In  what  follows,  a  dispersion 
analysis  is  performed  to  demonstrate  that  this  subtle  change  to  the  algorithm  results  in 
minimal  change  to  solutions  at  frequencies  that  are  well  resolved  by  the  grid  spacing, 
while  significantly  increasing  the  attenuation  at  higher  frequencies. 

A  simple  homogeneous  medium  with  uniform  density,  sound  speed  and  viscosity  is 

assumed.  In  this  case  the  analytic  solution  has  the  forme'c**r_fflt),  so  that  solutions  of  the 
form 


r  A;  «(£,.  /A  +  kjvA  -coNAt ) 

fj.,n  =  F0e 


(25) 


may  be  postulated  for  the  pressure  and  particle  velocities  in  each  direction.  The  subscripts 
indicate  the  grid  position  and  kr  and  kz  are  the  finite  difference  approximations  to  the 
exact  horizontal  and  vertical  wavenumbers.  The  numerical  approximation  to  the 
wavenumber  differs  from  the  analytic  value  of  k  by  an  amount  that  depends  on  the 
spatial  and  temporal  discretizations,  the  direction  of  propagation,  as  well  as  the  frequency 
and  sound  speeds.  If  0  is  the  propagation  angle  with  respect  to  the  radial  axis,  then 
kr  =  kcosO  and  k?  =  ksinO  For  simplicity,  the  dispersion  relation  is  derived  here  for  2D 
linear  propagation  in  Cartesian  coordinates,  i.e.  1/r  terms  are  neglected. 


As  in  de  Groot-Hedlin  et.al.  (2011),  the  following  notation  is  adopted 

St  =  sin(ojAf)/At  Sr  =  2sin(&rA/2)/A:  S.  =  2sin(/c.A /2)/ A-  (26) 

In  the  limit  of  small  Aand  At,  these  values  approach  to,  kr,  and  kz,  respectively  and 
(s;+s2z)=k\  Using  this  notation,  the  discretized  and  stabilized  linear  equations  in  2D 
are 


SP  =pc2(SV  +SV  ), 

to  /0  V  T  r  z  ZnJ’ 


(27) 
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(28) 


K5=^-  —{S2 + S2)\  ieiaAtV  +2 At 5^  | , 
r°  '  P0  P0  r  2  l  r°  Pa  ) 


V  S  =S^-/^(S2+S2)\ieioAtV  +2At^|.  (29) 

°  ‘  p0  p0  r  2  \  °  Po) 

Equations  (27)  through  (29)  are  combined  to  arrive  at  the  numerical  dispersion  relation 
for  the  stabilized  system  of  equations  in  a  homogeneous,  absorptive  medium: 


2  ^2  At 

c  — - 

Po 


{S2+S2f-{c2-l-^eioAt 

r  z  p 

•  o 


S)(S2+S2)+S2  =  0. 


(30) 


For  high  spatial  and  temporal  sampling  rates,  this  equation  becomes  a  fourth  order 
equation  for  the  complex-valued  wavenumber  k. 


c2^k* 

Po 


c2(  1 


1/.ICO  -2  2  r\ 

+  (0  =0 

Ptf 


(31) 


The  wavenumber  is  related  to  the  frequency  dependent  wave  speed  cp  and  attenuation  a 
through  k[cu)  —  cojc-\- This  expression  may  be  compared  with  the  dispersion 
relation  in  a  model  with  constant  viscosity 


k2  = 


co 


c(l-ipta>  /  pc2) 


(32) 


(de  Groot-Hedlin  et.  al,  2011),  which  is  a  second  order  equation  for  the  wavenumber. 
Equation  (32)  describes  the  dispersion  in  a  viscous  model,  which  is  a  real,  physical 
phenomenon  that  occurs  in  any  dissipative  medium.  In  contrast,  Eq.  (30)  describes 
dispersion  effects  that  result  both  from  the  smoothing  applied  in  stabilizing  the  FDTD 
algorithm  and  from  approximations  that  stem  from  discretizing  the  continuous  governing 
equations.  For  very  small  time  steps  and  low  kinematic  viscosities,  both  Eqs.  (30)  and 
(31)  reduce  to  Eq.  (32). 

The  frequency  dependent  sound  speeds  and  attenuation  were  computed  for  a  viscous 
medium  with  a  realistic  sound  speed  profile;  this  is  the  same  sound  speed  model  used  in 
the  third  example  in  Section  4.  The  ambient  atmospheric  density  profile  was  also 
computed  from  this  profile  under  the  assumption  that  the  atmosphere  is  in  hydrostatic 
equilibrium.  Figure  2  shows  the  density  profile,  as  well  as  the  frequency  dependent  sound 
speed  and  attenuation  profiles.  The  solid  lines  show  the  sound  speed  profiles  computed 
using  Eq.  (32)  for  the  physical  dispersion.  As  indicated,  the  speed  of  infrasound  is  very 
weakly  frequency  dependent  in  nature,  except  at  very  low  densities.  As  shown,  sound 
speeds  are  dispersive  within  the  thermosphere  at  frequencies  above  5  Hz,  with  higher 
velocities  at  higher  frequencies.  The  attenuation  varies  approximately  as  the  square  of  the 
propagation  frequency. 
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The  dashed  lines  show  the  sound  speeds  and  attenuations  computed  using  Eq  (31),  for  the 
stabilized  FDTD  equations,  using  a  value  of  At  =  0.037.  This  sampling  rate  is  also  used 
for  the  third  example  in  Section  4.  As  shown,  the  2-step  stabilization  procedure 
introduces  higher  attenuation  values  at  frequencies  that  are  poorly  resolved  by  the 
sampling  rate,  as  required,  but  also  changes  the  propagation  speeds.  The  results  show 
that,  at  this  sampling  rate,  the  stabilization  method  does  not  change  the  physics  of 
nonlinear  propagation  at  frequencies  at  low  frequencies  but  that  at  thermospheric 
altitudes,  the  dispersion  that  results  from  stabilization  significantly  alters  sound  speeds 
and  increases  attenuation  at  frequencies  above  5  Hz. 

Higher  temporal  sampling  rates,  i.e.  smaller  At,  would  affect  the  attenuation  and  sound 
speeds  at  higher  altitudes.  Thus,  at  least  to  altitudes  of  120  km,  the  stabilization  scheme 
causes  minimal  change  in  sound  speeds  and  attenuation  at  frequencies  that  are  well 
resolved  by  the  temporal  sampling  rate,  but  higher  frequencies  are  strongly  attenuated. 
The  attenuation  shown  by  the  dashed  lines  represent  the  minimum  values  at  each 
frequency;  attenuation  values  derived  from  the  numerical  dispersion  relation,  Eq.  (30), 
are  higher.  The  degree  to  which  they  are  higher  depends  on  the  order  of  accuracy  used  for 
the  finite  difference  operators.  In  this  paper,  the  FDTD  operators  are  accurate  to  second 
order. 


Attenuation  (nepers/km) 


Figure  2.  Profiles  of  a)  density,  b)  sound  speed,  and  c)  attenuation  as  a  function  of 
altitude.  The  solid  lines  show  the  frequency  dependent  sound  speeds  and  attenuation 
derived  from  the  wave  equation  for  a  viscous  medium.  The  dashed  lines  show  the  sound 
speeds  and  attenuation  values  for  the  stabilized  FDTD  algorithm  at  a  range  of 
frequencies,  computed  using  a  value  of  dt  =  0.037  s.  Higher  temporal  sampling  rates,  i.e. 
smaller  At,  would  lead  to  enhanced  attenuation  values  at  higher  altitudes. 
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C.  Source  Fields 


For  each  model,  the  solution  is  initialized  at  t  =0  by  setting  the  acoustic  velocities  v,  and 
v2  to  zero.  For  linear  propagation,  the  pressure  perturbation  is  initialized  as 


e  Wwf  (.05^^  j  ft  <  4^, 

0  R>4w 


(33) 


where  R  =  yr2  +[z-zgf  is  the  distance  from  the  source  center,  located  at  altitude  z0  and 

radius  r=0,  and  w  is  the  wavelength  of  the  acoustic  signal  at  the  source  altitude,  at  the 
highest  frequency  of  interest.  The  source  function  is  thus  a  truncated  Gaussian  pressure 
perturbation;  the  truncation  serves  to  initialize  the  perturbations  to  zero  outside  the  source 
region.  For  nonlinear  propagation,  the  solution  is  initialized  with  a  density  perturbation  of 
the  form 


/>SM=< 


Ae  cos(Rtt / Hw) 
0 


R<4w 

R>Aw 


(34) 


where  the  amplitude  A  controls  the  degree  of  nonlinearity.  The  perturbed  pressure 
corresponding  to  the  density  perturbation  is  initialized  using  Eq.  14  prior  to  the  time 
stepping  loop. 

D.  Computation  Windows 

The  acoustic  fields  computed  within  the  FDTD  method  are  very  sparse,  i.e.,  the  pressure 
perturbations  and  particle  velocity  fields  are  zero  over  most  of  the  computational  domain. 
A  method  that  takes  advantage  of  this  sparseness  by  performing  computations  only  over 
non-zero  field  variables  would  significantly  reduce  the  computational  load.  The  approach 
followed  here  takes  advantage  of  the  fact  that  for  long-range  infrasound  propagation,  the 
field  variables  are  non-zero  over  only  a  segment  of  the  propagation  path  at  any  given 
time.  This  allows  computations  to  be  performed  over  a  smaller  section  of  the  entire 
FDTD  model  domain  that  moves  along  with  the  pressure  pulses. 

An  example  is  shown  in  Figure  3  for  the  sound  speed  profile  shown  in  Figure  2,  for  a 
source  at  30  km  altitude.  The  pressure  solutions,  scaled  by  the  reciprocal  of  the 
atmospheric  densities,  are  shown  at  times  T=720  s,  990  s,  and  1260  s.  Ray  solutions  are 
superimposed  on  the  FDTD  solutions.  The  vertical  lines  show  the  edges  of  the  windows 
over  which  FDTD  computations  are  performed  at  each  time  step.  Results  for  this  model 
will  be  examined  further  in  section  4. 

The  selection  of  the  window  edges  over  which  the  computations  are  performed  is  guided 
in  part  by  the  ray  solution  for  the  given  atmospheric  model.  For  a  given  time,  the  right 
edge  of  the  computational  domain  is  set  to  a  value  greater  than  the  maximum  range 
attained  by  the  ray  solution.  Similarly,  the  left  edge  is  set  to  a  value  less  than  the 
minimum  range  attained  by  the  ray  solution.  In  practice,  the  window’s  right  starts  at  a 
minimum  of  500  nodes  width  and  expands  to  the  right  at  a  velocity  given  by  Cright,  where 
Cnght  is  greater  than  the  sound  speed  at  the  source  altitude.  After  a  time  delay,  TC|,  the  left 
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edge  of  the  computational  regime  moves  at  a  velocity  given  by  Qeft,  where  Qeft  is  less 
than  the  sound  speed  at  the  source  altitude.  The  time  delay  allows  for  the  acoustic  pulse 
to  propagate  upward  and  out  of  the  computational  regime  after  a  single  reflection  at  the 
rigid  boundary,  and  is  given  by  Td=  (ztop+z 0)/  Cieft.  In  general,  the  window  range  is 
approximately  equal  to  its  total  height  so  the  reduction  in  computation  time  depends  on 
the  model  geometry.  For  the  examples  shown  in  Section  4,  the  use  of  computation 
windows  reduces  computation  time  by  a  factor  of  about  4  to  5. 


Range  from  Source  (km) 


Figure  3.  FDTD  Pressure  fluctuation  snapshots,  with  superimposed  ray  solutions. 
Pressure  fluctuations  for  the  sound  speed  profile  shown  in  Figure  2,  with  a  source  at  30 
km,  for  the  times  shown  in  each  frame.  The  pressure  fluctuations  are  scaled  by  the 
reciprocal  of  the  square  root  of  the  atmospheric  densities.  Ray  solutions  are  shown  by 
darker  lines  for  stratospherically  ducted  rays,  and  lighter  lines  for  thermo  spherically 
ducted  rays.  The  left  and  right  edges  of  the  FDTD  computational  window  are  shown  in 
solid  vertical  lines. 

4.  RESULTS  AND  DISCUSSION 

In  this  section,  the  nonlinear  FDTD  infrasound  propagation  solver  is  applied  to  three 
environmental  models  to  examine  how  source  amplitudes  affect  signal  strengths  and  peak 
frequencies  for  receivers  at  various  distances  from  the  source.  The  first  two  examples 
feature  simple  models  with  source  amplitudes  that  range  from  nearly  linear,  with 
amplitudes  that  are  a  very  small  fraction  of  the  ambient  pressure,  to  strong  shocks,  with 
source  pressures  larger  than  the  ambient  values.  In  the  third  example,  a  realistic 
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atmospheric  profile  is  used  to  investigate  nonlinear  effects  in  long-range  infrasound 
propagation  for  sources  that  lead  to  weak  shocks,  where  the  acoustic  pressure 
perturbations  are  less  than  the  ambient  value,  but  not  negligible. 

A.  Model  1:  A  Uniform  Halfspace 

The  first  example  involves  the  prediction  of  nonlinear  acoustic  signals  that  propagate 
through  an  inviscid  uniform  halfspace,  terminated  by  a  rigid  surface  along  the  bottom. 

The  purpose  of  this  example  is  to  examine  the  changes  in  waveforms  and  spectra  that 
result  from  varying  the  source  amplitudes  for  the  simplest  possible  model.  The  ambient 
density  and  pressure  are  set  to  their  average  sea  level  values  and  the  ambient  temperature 
is  set  to  15°C,  giving  a  static  sound  speed  of  340  m/s.  The  source  is  centered  at  an  altitude 
of  150  m,  equal  to  the  source  height  used  for  many  nuclear  tests  conducted  as  part  of 
Operation  Plumbbob  at  the  Nevada  Test  site  in  1957.  Its  shape  is  given  by  Eq.  (34), 
where  w  =  32  m,  which  is  equal  to  the  wavelength  at  a  frequency  of  10.6  Hz,  given  the 
sound  speed  of  340  m/s.  Fifty-eight  virtual  receivers  are  located  at  pressure  nodes  in  the 
center  of  cells  along  the  rigid  surface,  at  ranges  from  just  below  the  source  to  5.7  km 
from  the  source,  in  increments  of  100  m. 

The  numerical  propagation  modeling  code  was  run  for  five  source  amplitudes  with 
maximum  density  perturbations  ranging  from  0.1%  to  80%  of  the  ambient  density;  from 
Eq.  (14)  this  yields  pressure  perturbations  from  .14%  to  128%  of  the  ambient  pressure 
fields.  The  spatial  and  temporal  grids  for  this  model  depend  on  source  size,  with  finer 
discretizations  for  the  larger  amplitude  sources  to  retain  accuracy  as  the  responses  shift 
from  lower  to  higher  frequencies.  For  the  smallest  sources,  with  pressure  perturbations  of 
0.14%  and  1.4%  of  the  ambient  field,  the  spatial  and  temporal  discretizations  are  given 
by  A=1.6  m  and  dt  =  1.3  ms.  Given  the  model  sound  speed  of  340  m/s,  this  spatial 
sampling  rate  gives  10  gridpoints  per  wavelength  at  a  frequency  of  21  Hz.  For  sources 
with  maximum  pressure  perturbations  of  29%  and  61%  of  the  ambient  field,  A  =1  m  and 
dt  =  0.84  ms;  for  the  largest  source  A=0.8  m  and  dt  =  0.67  ms. 

Figure  4  shows  the  infrasound  responses  as  a  function  of  time  for  virtual  receivers  at 
varying  ranges,  for  the  smallest  and  largest  sources  modeled.  Figure  4a  shows  that  the 
response  to  the  low  amplitude  source  is  nearly  linear,  with  pulses  that  scale  as  with  the 
inverse  of  the  distance  from  the  source,  as  expected  in  a  3D  medium  (Mathews  and 
Walker,  1970),  and  with  shapes  that  do  not  vary  with  range.  However,  Figure  4b  shows 
that  the  shape  of  the  pressure  pulses  changes  significantly  as  it  propagates  from  a  very 
large  amplitude  source;  the  wavefronts  steepen,  as  is  characteristic  of  nonlinear  acoustic 
propagation  (Pierce,  1981),  and  the  trailing  negative  excursion  becomes  drawn  out  and  is 
lower  in  amplitude  than  the  positive  perturbation. 
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Time  (s) 

Figure  4.  Pressure  waveforms  for  nonlinear  propagation  in  a  uniform  halfspace. 

Pressure  waveforms  for  virtual  receivers  at  200  m  range  increments  along  the  rigid 
surface,  with  a  source  at  150  m  altitude  for  top)  a  small  source  with  a  maximum  pressure 
perturbation  of  0.14%  of  the  ambient  pressure  and  b)  a  large  source  with  a  maximum 
pressure  perturbation  of  128%  of  the  ambient  pressure.  The  dashed  lines  indicate  the 
amplitude  loss  that  would  result  for  a  linear  source  in  a  3D  medium.  Note  the  waveform 
steepening  as  the  pulse  travels  away  from  the  larger  source. 


Changes  in  waveform  shape  are  accompanied  by  a  change  in  frequency  content.  The 
spectrograms  in  Figure  5  show  how  the  waveform  amplitude  spectra  evolve  as  a  function 
of  range  from  the  source,  for  sources  of  varying  amplitude.  For  the  lowest  source 
amplitude  shown,  the  spectrum  of  the  pulse  is  approximately  Gaussian  in  shape  with  a 
single  peak  frequency  at  2.6  Hz,  and  the  frequency  content  does  not  vary  with  increased 
range  from  the  source.  With  increased  source  size,  the  amplitude  spectrum  becomes 
flatter  overall,  as  energy  shifts  from  lower  to  higher  frequency  sidelobes,  and  the  higher 
amplitude  sources  generate  signals  with  larger  numbers  of  narrow  sidelobes. 
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Figure  5.  Spectra  for  waveforms  in  a  uniform  halfspace  with  varying  source  amplitudes. 
Amplitude  spectra  for  pressure  waveforms  propagated  through  an  inviscid  halfspace, 
showing  the  spectral  variation  with  range  for  sources  of  increasing  amplitude.  From  top 
to  bottom,  the  ratio  of  the  maximum  pressure  perturbations  to  the  ambient  pressure  for 
each  source  was  a)  0.01,  b)  0.29,  c)  0.61  and  d)  1.28.  The  source  pressures  ofb),c),  and 
d)  are  29.2  db,  35. 7  db,  and  42.1  greater  than  for  a).  The  colorbars  are  shifted 
accordingly. 


As  energy  is  redistributed  between  frequencies  during  nonlinear  propagation  from  the 
source,  the  dominant  frequency  shifts  to  lower  values.  Figure  6  shows  the  dominant 
frequency  and  transmission  loss  (TL)  as  a  function  of  range  for  the  four  sources 
examined  in  Figure  5.  Figure  6a  shows  that  there  is  no  change  in  the  dominant  frequency 
for  the  lowest  source  amplitude  modeled,  but  as  the  source  pressures  increase,  the  shift  to 
lower  peak  frequencies  becomes  more  pronounced.  Figure  6b  indicates  that,  at  all  source 
sizes,  the  amplitude  varies  with  the  inverse  of  the  distance.  These  results  show  that 
nonlinearity  has  a  greater  effect  on  the  waveform  spectra  than  on  amplitudes  in  an 
inviscid  halfspace  model.  Once  the  pulse  pressure  amplitudes  have  dropped  to  the  linear 
regime,  the  form  of  the  frequency  spectrum  remains  unchanged  as  the  pulse  propagates 
further.  That  is,  the  spectra  retain  a  series  of  sidelobes  that  are  characteristic  of  their 
nonlinear  origins.  This  is  in  agreement  with  the  results  of  Castor  et  al.  (2004),  who  found 
that  the  spectra  of  hydroacoustic  waves  also  retain  the  nonlinear  signature  of  a  high 
amplitude  source  even  after  propagating  great  distances. 
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Range  from  Source  (km) 


Figure  6.  Variation  in  waveform  amplitude  and  peak  frequency  in  a  uniform  halfspace. 
Top)  The  peak  frequency  as  a  function  of  range  for  the  four  sources  examined  in  Figure 
5.  The  small  dots  show  the  dominant  frequency  for  the  smallest  source  amplitude.  Results 
for  the  larger  sources  are  shown,  by  increasing  source  amplitude,  by  triangles,  circles, 
and  squares,  respectively,  bottom)  The  transmission  loss  (TL)  as  a  function  of  range  for 
each  source  amplitude. 


B.  Model  2:  A  Bilinear  Halfspace 

The  second  example  examines  nonlinear  infrasound  propagation  into  a  shadow  zone  for  a 
model  with  a  bilinear  sound  speed  profile,  which  has  a  sound  speed  profile  given  by 

c(z )  =  c(0)  /  •> J\  +  2z7r  ,  where  R  is  the  radius  of  curvature  of  a  ray  at  the  source  altitude. 
An  asymptotic  solution  for  pressure  exists  for  this  type  of  sound  speed  profile  (Pierce, 
1981),  which  predicts  the  existence  of  a  creeping  wave  along  the  surface  within  the 
shadow  zone.  This  creeping  wave  has  an  attenuation  term  that  is  proportional  to  the  third 
root  of  the  frequency. 

The  sound  speed  profile  for  this  model  is  given  by  c(z)=  340m  / s  x  a/i  +  z  /  39000 , 

which  gives  a  nearly  linear  sound  speed  decrease  of  4.2  m/s  per  kilometer  altitude.  This  is 
similar  to  the  rate  of  decrease  resulting  from  the  environmental  lapse  rate  of  6.5°  C  per 
km,  the  rate  of  temperature  decrease  with  altitude  in  a  stationary  atmosphere  from  sea 
level  to  1 1  km.  The  source  is  identical  to  that  of  the  previous  model,  and  the  ambient 
densities  and  pressures  are  set  to  their  sea  level  values.  The  nonlinear  code  was  run  for 
four  source  amplitudes  with  maximum  density  perturbations  ranging  from  1%  to  80%  of 
the  ambient  density,  equivalent  to  maximum  pressure  perturbations  from  1.4%  to  128% 
of  the  ambient  pressure  fields.  Again,  viscosity  is  neglected  as  it  is  negligible  for  the 
densities  and  frequencies  considered  here. 

Figure  7  shows  the  sound  speed  profile  and  ray  solution  for  this  model,  along  with  the 
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pressure  field  computed  for  the  smallest  source  at  time  60  s  after  the  shot.  The  ray 
solution  shows  a  shadow  zone  starting  at  a  range  of  approximately  5  km  from  the  source, 
but  the  FDTD  solution  indicates  that  acoustic  energy  penetrates  well  into  the  shadow 
zone.  The  nonlinear  FDTD  solution  was  run  with  80  receivers  along  the  rigid  surface  at 
intervals  of  250  m,  beginning  just  below  the  source. 


Sound  Speed  (m/s) 


Range  from  Source  (km) 


Figure  7.  a)  sound  speed  profile  for  a  bilinear  halfspace,  b)  Ray  solution  with 
superimposed  FDTD  pressure  solution.  The  bilinear  sound  speed  profile  used  for 
Example  2  is  shown  at  left  and  the  ray  solution  is  at  right,  superimposed  on  the  FDTD 
solution  computed  at  60  seconds.  The  circle  at  150  m  altitude  at  left  shows  the  source 
location.  Triangles  mark  the  locations  of  the  virtual  receivers.  The  ray  solution  indicates 
a  shadow  zone  beyond  approximately  7  km  from  the  source,  but  the  FDTD  solution 
indicates  that  infrasound  penetrates  into  this  region. 


The  spectrograms  in  Figure  8  show  the  evolution  of  the  amplitude  spectra  with  increasing 
range  for  this  sound  speed  profile  for  sources  of  varying  size.  Compare  this  to  Figure  5, 
which  shows  equivalent  results  for  a  halfspace  model.  As  for  the  halfspace,  the  spectra 
become  progressively  more  scalloped  as  the  source  amplitude  increases,  with  a  shift  to 
higher  frequency  sidelobes.  Figure  8  shows  that  the  higher  frequencies  are  more  severely 
attenuated  within  the  shadow  zone  than  low  frequencies,  in  agreement  with  the 
asymptotic  solution  for  a  bilinear  model. 

Figure  9  shows  the  dominant  frequency  and  TL  as  a  function  of  range  for  these  sources. 
The  dominant  frequency  shifts  to  lower  values  for  propagation  into  a  shadow  zone  at  all 
the  source  amplitudes,  as  indicated  in  Figure  9a.  In  this  case,  the  peak  frequencies 
continue  to  decrease  with  range  even  after  the  pressure  fluctuations  have  decreased  to  the 
linear  regime,  as  may  be  expected  for  propagation  of  a  creeping  wave  within  a  shadow 
zone.  Figure  9b  shows  that  the  source  amplitude  does  not  have  a  significant  effect  on  the 
TL,  even  within  the  shadow  zone.  However,  the  sound  speed  profile  has  a  strong  effect 
on  the  TL  at  all  source  amplitudes,  with  greater  losses  in  signal  strength  within  the 
shadow  zone  than  for  a  uniform  halfspace  with  an  identical  source. 
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Figure  8.  Spectra  for  waveforms  in  a  bilinear  halfspace  with  varying  source  amplitudes. 
As  for  Figure  5,  but  for  a  model  with  the  sound  speed  profile  shown  in  Figure  7. 


Range  from  Source  (km) 


Figure  9.  Variation  in  waveform  amplitude  and  peak  frequency  in  a  bilinear  halfspace. 
As  for  Figure  6,  but  for  a  model  with  the  sound  speed  profile  shown  in  Figure  7.  The 
dashed  line  in  the  lower  panel  indicates  the  TL  expected  for  a  uniform  halfspace. 
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C.  Model  3:  Propagation  through  the  stratosphere  and  thermosphere 

In  the  final  example,  infrasound  propagation  is  computed  for  a  viscous  medium  with  a 
realistic  static  sound  speed  profile,  to  examine  the  effects  of  nonlinearity  on  propagation 
within  the  thermosphere,  where  the  ambient  density  is  very  low,  or  at  caustics  within 
stratospheric  ducts.  Both  the  linear  and  nonlinear  codes  are  applied  to  this  model.  The 
sound  speed  and  density  profiles  used  for  this  model  are  shown  in  Figure  2.  Virtual 
infrasound  receivers  were  placed  near  the  rigid  surface  at  intervals  of  4  km  to  a  range  of 
452  km. 

The  source  altitude  was  set  to  30  km,  approximately  the  height  of  the  Chelyabinsk  bolide 
(de  Groot-Hedlin  and  Hedlin,  2014)  and  also  that  of  some  shots  used  for  a  high-altitude 
infrasound  calibration  experiment  (Herrin  et.al,  2008).  The  source  shape  is  given  by  Eq. 
(34)  with  w  =  606  m.  For  the  nonlinear  sources,  the  maximum  density  perturbations  are 
set  to  1%  and  4%  of  the  maximum  density  perturbations  at  the  source  altitude,  equal  to 
1.4%  and  5.6%  of  the  ambient  pressure  fields  at  30  km.  Although  these  sources  have 
much  lower  amplitude  than  in  the  previous  examples,  the  total  spatial  volume  is  about 
6800  times  greater,  leading  to  a  much  lower  peak  frequency.  The  discretization  levels 
were  set  to  A  =52  m  and  dt  =  0.037  s  for  linear  propagation  modeling,  and  to  A  =  26  m 
and  dt  =0.018  s  for  both  nonlinear  propagation  tests. 

Akintunde  and  Petculescu,  (2014)  suggest  that  it  may  be  inaccurate  to  treat  the  upper 
altitude  as  a  continuum  fluid  at  high  frequencies;  their  results  suggests  the  use  of  the 
Navier-Stokes  equation  to  describe  infrasound  propagation  may  be  inaccurate  at  altitudes 
above  120  km.  and  frequencies  above  1  Hz.  Accordingly,  the  model’s  maximum  altitude 
was  set  to  120  km.  The  ray  solution  and  the  FDTD  pressure  solutions  are  shown  in  Figure 
3  at  several  time  points.  The  maximum  sound  speed  within  the  stratosphere  is  less  than 
the  sound  speed  ground  level,  so  the  ray  solution  predicts  that  infrasound  would  be 
trapped  within  a  stratospheric  duct  and  not  recorded  at  the  surface.  It  predicts  to  a  shadow 
zone  from  approximately  100  km  to  nearly  300  km  from  the  source,  where  the  first 
thermospheric  returns  arrive.  The  ray  solution  also  predicts  the  presence  of  caustics 
within  the  stratosphere  at  ranges  of  approximately  125  km  and  250  km  from  the  source. 
The  FDTD  solution  confirms  the  presence  of  caustics  within  the  stratosphere. 

Figure  10  shows  pressure  responses  synthesized  using  the  nonlinear  algorithm  for  the 
smaller  source,  scaled  by  the  distance  from  the  source.  The  ray  arrival  times,  indicating 
only  arrivals  that  reach  the  ground,  are  shown  for  comparison.  As  shown,  the  arrival 
times  for  the  FDTD  waveforms  agree  with  the  ray  travel  times  for  direct  and 
thermospheric  arrivals,  however  the  FDTD  waveforms  also  include  low  amplitude 
stratospheric  arrivals.  Comparing  this  to  Figure  3,  it  can  be  seen  that  these  stratospheric 
arrivals  are  strongest  where  the  rays  nearly  reach  the  ground. 

The  spectrograms  for  the  linear  and  nonlinear  solutions  are  compared  in  Figure  11.  The 
colorscale  of  the  linear  modeling  result  is  arbitrary,  as  doubling  the  source  size  leads  to 
waveforms  with  exactly  twice  the  amplitude  with  no  change  in  the  waveform  shape.  The 
spectral  amplitudes  for  the  linear  results  were  therefore  uniformly  shifted  to  agree  with 
the  results  for  the  smaller  nonlinear  source  near  the  source  origin.  As  shown,  there  is  little 
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difference  in  the  frequency  content  for  the  direct  arrivals,  at  ranges  from  0  to  just  over 
100  km.  A  presence  of  a  single,  small  sidelobe  in  the  spectrogram  for  the  largest 
amplitude  source  indicates  a  weakly  nonlinear  response  in  the  near-source  region. 


Reduced  Time  (s)  (V  =  350  m/s) 


Figure  10.  Waveforms  computed  by  FDTD,  with  predicted  ray  arrival  times.  Waveforms 
synthesized  by  the  nonlinear  FDTD  code,  scaled  by  the  source  distance,  are  shown  as  a 
function  of  source  range.  The  waveforms  were  synthesized  for  the  sound  speed  and 
density  profdes  shown  in  Figure  2.  The  arrival  times  predicted  by  rays  are  shown  by 
dots,  showing  only  arrivals  that  reach  the  ground.  Small  amplitude  stratospheric  arrivals 
can  be  seen  within  the  shadow  zone  predicted  by  the  ray  solution.  Spectrograms  in 
Figure  11  are  computed  over  the  time  windows  shown  in  gray. 


The  low  frequency  stratospheric  arrivals  within  the  ‘shadow  zone’  between  100  and  280 
km  have  much  lower  frequency  content  that  the  direct  arrivals  and  are  nearly  identical  for 
the  linear  and  nonlinear  computations.  The  low  frequency  content  is  characteristic  of 
creeping  waves  and  does  not  result  from  the  shear  viscosity  of  the  medium.  To  verify 
this,  the  linear  code  was  rerun  for  the  same  sound  speed  and  density  profiles,  but  for  an 
inviscid  medium.  The  results  confirm  that  shear  viscosity  has  no  effect  on  the 
stratospheric  arrivals.  The  similarity  of  the  spectral  results  at  this  range  indicates  that  the 
existence  of  caustics  in  the  solution  does  not  cause  significant  changes  to  the  signals  for 
relatively  low  amplitude  sources. 

Thermospheric  arrivals  that  appear  at  ranges  beyond  250  km  have  significantly  lower 
frequency  content  than  the  direct  arrivals.  Comparisons  between  the  viscous  and  inviscid 
linear  solutions  indicate  that  viscosity  significantly  attenuates  thermospheric  returns  at 
frequencies  over  0.2  Hz.  Comparisons  between  the  linear  and  nonlinear  solutions  show 
clear  differences,  more  so  than  for  the  direct  arrivals.  For  infrasound  propagation 
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through  the  thermosphere,  computed  for  finite  amplitude  sources,  there  is  an  opposing 
effect  between  atmospheric  viscosity,  which  lowers  the  frequency  of  the  arrivals  and 
nonlinear  effects,  which  steepen  the  wavefronts  of  the  signals,  shifting  energy  to  higher 
frequency  sidelobes.  The  spectrograms  of  the  thermospheric  arrivals  show  a  shift  of 
acoustic  energy  to  higher  frequency  sidelobes,  indicating  a  relatively  high  degree  of 
nonlinearity  even  for  the  low  amplitude  source.  The  maximum  amplitude  for  the 
thermospheric  arrivals  from  the  smaller  source  is  about  0. 1  Pa,  for  the  larger  source  it  is 
0.5  Pa,  suggesting  that  any  detectable  thermospheric  returns  are  altered  by  nonlinear 
effects. 
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Figure  11.  Spectra  for  waveforms  in  a  realistic  atmospheric  model.  Amplitude  spectra 
for  propagation  through  the  sound  speed  profile  shown  in  Figure  3,  showing  the  spectral 
variation  with  range  for  sources  of  increasing  amplitude.  The  top  panel  presents  the 
results  for  linear  propagation  modeling,  the  lower  two  present  the  results  for  nonlinear 
propagation.  For  the  middle  and  bottom  panels,  the  ratio  of  the  maximum  pressure 
perturbations  to  the  ambient  pressure  at  30  km  was  middle)  0.014,  bottom)  0.056.  The 
source  pressure  of  the  larger  source  is  12.1  db  larger  than  for  the  smaller  source.  The 
colorbar  is  shifted  accordingly.  The  amplitudes  for  the  linear  results  are  arbitrary  and 
were  scaled  to  agree  with  the  results  for  the  smaller  source  near  the  source  origin. 


Figure  12  shows  the  dominant  frequency  and  TL  as  a  function  of  range  for  this  model,  for 
the  linear  and  nonlinear  code.  There  are  insignificant  differences  in  the  peak  frequencies 
for  the  linear  and  nonlinear  results  up  to  200  km  range.  At  greater  ranges,  the  peak 
frequencies  for  the  smaller  source  are  similar  to  those  for  the  linear  model;  the  peak 
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frequencies  are  lower  for  the  larger  source.  As  for  the  simpler  models,  the  sound  speed 
profile  has  a  significant  effect  on  the  TL,  but  the  relative  source  amplitude  does  not. 


Figure  12.  Peak  frequencies  and  TL  as  a  function  of  range  for  a  realistic  atmosphere. 
Top)  The  peak  frequency  as  a  function  of  range  for  the  cases  examined  in  Figure  11.  The 
solid  line  shows  the  dominant  frequency  for  the  linear  modeling  results.  Results  for  the 
nonlinear  sources  are  shown  by  circles  for  the  smaller  source  and  squares  for  the  larger 
source,  bottom)  The  transmission  loss  (TL)  as  a  function  of  range  for  each  source 
amplitude. 

5.  CONCLUSIONS 

FDTD  algorithms  to  compute  inffasound  propagation  through  a  viscous  medium  can  be 
unstable  at  very  low  atmospheric  density.  A  computationally  efficient  method  to  stabilize 
FDTD  codes  has  been  developed  for  both  linear  and  nonlinear  infrasound  propagation 
modeling.  The  stabilization  is  accomplished  by  breaking  up  the  computation  of  the 
particle  velocities  into  two  steps:  the  acoustic  term  is  added  in  the  first  step  and  the 
viscosity  operator  acts  upon  the  resulting  term.  This  scheme  results  in  a  subtle  alteration 
of  the  governing  equations,  and  thus  alters  the  dispersion  characteristics.  This  dispersion 
is  quantifiable  and  it  is  shown  that  sound  speeds  and  attenuation  are  unchanged  at 
frequencies  that  are  well  resolved  by  the  temporal  sampling  rate,  but  higher  frequencies 
are  strongly  attenuated,  thus  suppressing  instability. 

Applying  the  FDTD  algorithms  to  strong  shocks  in  inviscid  models  and  weak  shocks  in 
viscous  models  shows  that  the  source  amplitudes  have  a  stronger  effect  on  the  dominant 
frequency  than  on  the  TL.  However,  the  amplitude  and  spectral  content  of  infrasound 
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signals  are  largely  controlled  by  the  sound  speed  profile,  indicating  that  the  empirical 
source  yield  relations  used  for  nuclear  monitoring  that  are  based  on  either  the  signal 
amplitude  or  peak  frequency  must  be  amended  for  the  effective  sound  speed  profiles.  The 
fact  that  TL  is  not  strongly  dependent  on  source  amplitude  suggests  that  the  expressions 
for  TL  developed  by  Le  Pichon  et.al.  (2012)  for  a  wide  range  of  atmospheric  conditions 
can  be  applied  to  propagation  of  weak  shocks  too.  More  work  is  required  to  examine  the 
TL  for  infrasound  propagation  through  the  thermosphere  for  sources  on  the  same  order  of 
magnitude  as  the  ambient  pressure. 

The  results  from  applying  the  stabilized  FDTD  code  to  infrasound  propagation  through  a 
realistic  sound  speed  and  density  profile  show  that  thermospheric  arrivals  are  strongly 
affected  by  nonlinear  effects,  more  so  than  for  arrivals  near  the  source.  Energy  is  shifted 
from  lower  to  higher  frequencies  for  thermospheric  returns,  even  for  signal  amplitudes  of 
less  than  1  Pa.  This  indicates  that  nonlinear  effects  are  significant  for  any  detectable 
thermospheric  returns. 

A  possible  extension  of  this  work  is  to  apply  the  stabilized  nonlinear  FDTD  code  to 
improve  our  understanding  of  acoustic  propagation  in  the  thermosphere.  Although  source 
yield  estimates  can  be  made  for  infrasound  returns  from  the  thermosphere,  and  indeed, 
thermospheric  returns  are  often  the  only  signals  observable  at  distance  recording  sites  (Le 
Pichon  et  al.,  2013;  de  Groot-Hedlin  and  Hedlin,  2014),  infrasound  propagation  in  the 
thermosphere  remains  poorly  understood.  Acoustic  absorption  and  dispersion  are  poorly 
constrained  at  altitudes  above  90  km  (Sutherland  and  Bass,  2004)  and  the  treatment  of  the 
upper  atmosphere  as  a  continuum  fluid  has  been  called  into  question  (Akintunde  and 
Petculescu,  2014).  Comparisons  between  observations  and  the  results  of  waveform 
modeling  may  resolve  some  of  these  unknowns. 
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